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ABSTRACT 

A  method  of  digital  computer  oriented  analysis  of  a 
Mason  flow  graph  is  presented.   Using  this  method  in  connec- 
tion with  a  computer  program  implementing  the  algebraic 
manipulation  language  FORMAC  can  be  a  tool  in  expanding  large 
determinants  with  symbolic  polynomial  elements  and  thereby 
finding  the  characteristic  equation  and  transfer  functions  of 
a  linear  time-invariant  system  in  symbolic  form. 

Also  a  new  approach  to  the  root  locus  method  is  shown, 
using  a  FORMAC  program.   The  advantages  over  the  conventional 
root  locus  method  are  discussed. 

Areas  of  possible  future  use  of  FORMAC  in  algebraic 
problems  of  control  theory  are  discussed. 
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I.   INTRODUCTION 

A.  THE  NATURE  OF  THE  PROBLEM 

Consider  a  linear  time-invariant  system  whose  parameters 
are  known  only  symbolically.   In  the  analysis  of  the  system 
often  it  becomes  necessary  to  derive  certain  system  functions 
which  are  polynomials  in  the  Laplace  variable  s. 

The  "manual"  derivation  of  these  functions,  even  for 
small  systems,  becomes  a  task  which  is  not  only  time  consum- 
ing but  unattractive,  and  although  conceptually  simple,  is 
also  a  source  of  errors  which  are  difficult  to  detect. 

The  enormous  work  involved  in  the  deriving  of  these 
functions  for  medium  or  large  systems  can  be  the  reason  for 
not  attempting  a  solution  to  the  problem  in  general  algebraic 
form. 

There  are  also  possible  other  algebraic  problems  where 
the  "manual"  derivation  is  tedious. 

The  problem,  therefore,  is  to  devise  methods  of  using 
the  digital  computer  to  perform  those  algebraic  manipulations. 

B.  METHOD  OF  SOLUTION 

The  first  specific  problem  to  be  considered  will  be  the 
evaluation  of  the  characteristic  equation  of  a  linear  system. 

To  solve  this  problem  it  will  be  necessary  to  evaluate 
the  determinant  of  a  square  matrix  whose  elements  are  generally 
polynomials  in  the  variable  s,  derived  from  the  differential 
equations  of  the  system  by  Laplace  transform. 

Two  different  methods  of  evaluating  a  determinant  are 
applied  here.   The  first  method  is  based  on  the  theory  of 
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graphs,  specifically  on  the  properties  of  a  Mason  flow  graph 
(1) .   It  will  be  assumed  that  the  reader  is  familiar  with  the 
Mason  gain  formula. 

Other  graph  methods  are  those  of  Coates  and  the  modified 
method  of  Chan  and  Bapna  [Ref .  2] .   These  methods  are  not 
based  on  a  Mason  flow  graph.   A  flow  graph   analysis  using  a 
Mason  flow  graph  and  based  on  a  topological  development  is 
given  by  Dunn  and  Chan  [Ref.  3]. 

The  flow  graph  analysis  used  here  is  not  pure  topological, 
because  the  loops  and  paths  are  not  found  from  the  topological 
properties  of  the  graph  but  by  search  methods. 

A  comparison  of  flow  graph  techniques  and  another 
different  method  by  Chan  and  Mai  are  found  in  [Ref.  2]. 

The  basic  feature  of  the  mechanized  graph  analysis  used 
here  is  the  process  of  identifying  closed  unidirectional  paths 
called  "loops"  among  the  connections  of  the  graph.   Then  the 
set  of  non-touching  loops  is  found  which  is  then  used  to 
evaluate  the  determinant. 

The  second  method  used  here  for  the  evaluation  of  a 
determinant  is  using  a  fraction-free  algorithm.   Lipson  [Ref. 
4]  uses  this  algorithm  in  the  evaluation  of  a  set  of  linear 
algebraic  equations. 

The  basic  problem  in  the  evaluation  of  the  determinant 
Det (A)  is  the  identification  of  those  elements  in  the  square 
matrix  A  whose  products  form  a  term  in  the  determinant.   There 
are  n!  terms  in  the  determinant  of  a  square  matrix  of  order 


n;  most  of  them  are  often  equal  to  zero  in  a  linear  system  of 
the  kind  considered  in  control  theory. 

The  computer  program  DTERM  (Program  1)  shows   in  its 
algebraic  part  how  products  of  polynomials  can  be  processed 
without  using  excessive  space  (which  was  the  main  limitation 
imposed) . 

The  first,  non-algebraic  part  of  DTERM  uses  the  basic 
definition  of  a  determinant  to  find  the  non-zero  products  in 
the  determinant  by  forming  all  possible  permutations  of  the 
column  indices  of  the  elements  of  the  matrix  A.   This  method 
was  used  only  to  simulate  the  presence  of  a  result  of  the 
flow  graph  analysis,  which  was  not  programmed. 

For  matrices  of  order  5  it  was  found  that  the  execution 
time  for  investigation  of  all  120  possible  products  of  5 
elements  was  6  see,  for  order  8  the  execution  time  was  27  sec 
and  for  order  9  the  execution  time  was  3  min  43  sec.   These 
times  include  also  other  statements. 

The  fraction-free  algorithm  is  given  here  as  a  different 
possible  method  to  evaluate  the  determinant  of  a  matrix,  using 
only  single  letters  as  symbolic  elements  followed  by  a  gradual 
replacement  by  the  longer  polynomials.   The  latter  process  is 
the  problematic  one  because  of  the  enormous  space  requirement 
which  can  be  generated.   It  would  probably  be  similar  to  the 
method  used  in  the  algebraic  part  of  DTERM. 

For  the  other  algebraic  problems  mentioned  here  short 
example  programs  in  FORMAC,  an  algebraic  manipulation  language, 


are  given.   These  programs  show  the  possible  use  of  FORMAC  in 
the  solution  of  those  problems;  they  are  not  optimal. 

C.    THE  ALGEBRAIC  MANIPULATION  LANGUAGE  "FORMAC" 

The  algebraic  manipulation  must  be  performed  with  the 
aid  of  a  special  tool:   an  algebraic  manipulation  language. 
Here  "FORMAC"  was  used  for  this  purpose. 

FORMAC  stands  for  FORmula  MAnipulation  Compiler.   It  was 
developed  by  the  IBM  Corporation  as  "an  experimental  system  to 
provide  engineers  with  the  capability  of  handling  formal 
mathematical  expressions  on  a  computer " [5] . 

Examples  of  application  of  FORMAC  are  given  in^"^  and  '  ', 
Reference  to  bibliographies  is  made  in  ^°' . 

The  most  recent  form  of  FORMAC  is  available  in  the  PLl- 
Version.   The  language  PL/1  is  a  proper  subset  of  FORMAC,  it 

is  therefore  possible  to  use  the  entire  PL/1  capability  toget- 

(9) 
her  with  FORMAC.   The  basic  reference  manual  for  PL/1  v 

( 6 ) 
together  with  the  respective  manual  for  FORMAC  l    at  this 

time  (1970)  is  all  the  user  must  be  referred  to.   The  average 
FORTRAN  user  can  however  in  a  few  days  acquire  effective  know- 
ledge of  both  PL/1  and  FORMAC  in  order  to  communicate  with  the 
FORMAC  language. 

FORMAC  enables  the  user  to  analyze  and  synthesize  general 
algebraic  expressions. 

The  work  with  FORMAC  proceeds  basically  in  such  a  way 
that  the  user  writes  a  program  in  the  FORMAC  language,  follow- 
ing the  route  of  analysis  of  normal  "manual"  algebraic 

manipulation . 
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The  capabilities  of  FORMAC  can  be  divided  into  the 
following  groups: 

-  analysis:    test  for  mathematical  identity,  finding 
the  coefficient  of  an  expression  inside 
a  parent  expression,  finding  the  high- 
est or  lowest  power  of  one  expression 
inside  a  parent  expression,  finding 
the  numerator  or  denominator  of  a 
rational  expression,  finding  the  main 
operator  of  an  expression  (+,-,:,x, 
and  so  on) ,  finding  the  number  of  terms 
of  an  expression  (factors  or  summands) , 
selection  of  a  specified  term  of  an 
expression. 

-synthesis:    assignment  of  an  algebraic  expression 
to  a  FORMAC  variable,  user  defined 
function  processing,  conversion  routines, 
well  developed  transit  between  "normal" 
values  and  PL/1  variables  and  FORMAC 
variables,  multiplication,  division, 
addition  and  subtraction,  differentia- 
tion with  respect  to  specified  and 
unspecified  arguments,  use  of  trigono- 
metric and  natural  functions  with 
numeric  and  symbolic  arguments,  combin- 
atorial and  factorial  integer  functions. 
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■support:      simple  printout,  editing,  control  over 
level  of  expansion  and  evaluation. 
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II.   METHOD  OF  COMPUTER  ANALYSIS 
FOR  A  MASON  FLOW  GRAPH 

In  this  chapter  a  procedure  will  be  described  which  can 
be  implemented  using  a  digital  computer.   The  basic  operations 
performed  by  the  computer  are  list  searches  and  logic  oper- 
ations on  binary  matrices. 

The  method  produces : 

-  lists  of  node  numbers  which  form  a  specific  loop. 

-  lists  of  associated  branch  gains,  so  loop  gains 
can  be  evaluated. 

-  a  list  containing  the  names  of  touching  loops,  so 
sets  of  touching  loops  can  be  identified. 

-  lists  of  node  numbers  which  form  a  path  from  a 
specified  node  to  another  specified  node.   All  the 
lists  identify  all  the  paths  present  between  these 
nodes . 

-  lists  of  associated  branch  gains,  so  the  different 
path  gains  can  be  evaluated. 

-  a  list  of  branch  point  nodes  which  identify  the 
points  where  the  different  paths  branch. 

First  the  terms  used  will  be  defined,  then  the  method  will 
be  described,  using  examples.  The  computer  implementation  will 
be  discussed. 

A.    DEFINITIONS 

1.   The  Mason  Flow  Graph 

Here  the  definitions  and  notation  of  Ward  and  Strum [1] 
will  be  used. 
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The  Mason  flow  graph  consists  of  a  set  of  nodes, 
which  are  here  arbitrarily  numbered,  and  a  set  of  branches, 
which  have  a  direction  and  a  weight  associated  with  them.   The 
branches  connect  the  nodes  in  such  a  way  that  at  least  every 
node  is  connected  with  one  other  node;  self  loops  are  not 
present. 

2*   Input-  and  Output-Nodes,  Branchpoint  Node 

An  input  node  has  no  incoming  branches.   Here  an  out- 
put node  has  one  incoming  branch.   A  branchpoint  node  is  a  node 
with  at  least  two  outgoing  branches.   From  a  branchpoint  node 
there  are  at  least  two  paths  leading  to  a  common  end  node. 

End-  and  starting  nodes  are  the  end-  and  the  starting 
point  respectively  of  paths. 
3.   Input  List  LIST1 

The  input  data  are  prepared  in  the  form  of  an  input 
list,  which  is  called  "LIST1".  Each  item  in  LIST1  has  three 
parts:   two  node  numbers  and  a  branch  designator. 

These  definitions  are  illustrated  in  Fig.  1  and  2. 
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Fig.  1 
Nodes  3,  2,  6  are  input  nodes.   This  graph  contains  no 
branchpoint  nodes.   The  other  nodes  are  normal  nodes, 
except  node  9,  which  is  output  node. 


? f 


Fig.  2 
In  this  graph  node  4  is  (also)  a  branchpoint  node. 
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LIST1 
Column  \  Column  2  Column  3 

node  number  k      branch  designator  A    node  number  1 

from  via  to 

elements  in  Col.  1  are  referred  to  as  LIST1  (l,k) 
elements  in  Col.  3  are  referred  to  as  LIST1  (3,k) 

LIST1  constitutes  a  complete  listing  of  the  graph. 
Columns  1  and  3  contain  numbers,  Column  2  contains  an  alpha- 
numeric designator  in  the  form  "letter-number",  "letter"  or 
"letter-letter". 

The  order  of  the  listing  is  arbitrary;  it  is  recommend- 
ed to  bring  the  list  in  at  least  approximate  numerical  order  in 
Column  1  for  improved  error  detection. 
4.   LIST2  and  LIST3 

LIST2  is  identical  to  LIST1  except  Columns  1  and  3 
are  exchanged  and  then  Column  1  in  LIST2  is  arranged  in 
numerical  order.   LIST3  has  the  same  structure  as  LIST1,  but 
the  contents  have  changed. 

5-   Topological  Matrices 
a.   Loop  Matrix 

A  loop  matrix  L  is  defined.   The  elements  of  L 


are  ajk. 


ajk  "  1       f°r  node  k  being  in  loop  j 
=  0       o.w. 

j=l,2,3, . . ,nl,  nl  =  number  of  loops 
k-1,2,3, . . ,nn,  nn  =  number  of  nodes 
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LT   = 


110    0    0  0    0 

0    0    110  0    0 

0    0    0    1    1  1    1 

0    0    0    0    0  11 


>T 


1 

1 

1 

1 

1 

1 

1 

0 

0 

1 

0 

0 

1 

0 

0 

1 

1 

0 

1 

1 

1 

Fig.    3 
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b.   Path  Matrix  P 

A  path  matrix  P  is  defined.   The  elements  of  P 
are  bmk- 

mk  =  -1       for  node  k  being  in  path  m 
=  °       o.w. 

m  =  l,2,3,...nn,  nn  =  number  of  nodes 
c   Touching  Loops  Matrix  LT 

A  matrix  is  defined  which  contains  information 
about  loops  which  touch  in  the  Mason  flow  graph.   The  elements 
of  this  matrix  LT  are  dn  . . 

dkj  =  1       for  loop  k  touches  loop  j;  that  is  loops 
k  and  j  have  at  least  one  node  in 
common . 
-  0       o.w. 

k,  j        l,2,3,...nl,  nl  =  number  of  loops 
d.   Path  -  Loop  Touching  Matrix  PL 

A  matrix  PL  is  defined  which  indicates  which        I 

paths  touch  which  loops.   The  elements  of  PL  are  e 

km* 

ekm  =  -1       for  looP  k  touches  path  m 
=  0       o.w. 

k=l,2,3, . . .nl,  nl  =  number  of  loops 
m=l,2,3, . . ,np,  np  =  number  of  paths 
e.   Control  Vectors 

Certain  control  vectors  are  defined  to  check  on 
specified  conditions. 
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(1)  Loop  Search  Control  Vector  LC 
The  element  of  LC  are  c-.^. 

c,,  =  1       for  every  input  and  output  node 

=  1       for  every  other  node  used  at  least  once 

in  the  loop  search 
=  0       o.w. 

k=l,2 , 3 , . . .nn,  nn  =  number  of  nodes 

(2)  Multiple  Node  Control  Vector  LIST1(4) 
LISTl(4,k)  =  0   for  node  LISTl(l,k)  being  a  single 

node  or  the  mth  node  in  a  series  of 
m  equal  node  numbers 
=  1   for  the  first  m  -  1  node  numbers  in 
a  sequence  of  m  equal  node  numbers 
k=l,2,3/ . . .nn,  nn  =  number  of  nodes 

(3)  Branchpoint  Control  Vector  BC 

This  vector  contains  information  about  nodes 
which  are  branchpoint  nodes.   The  information  is  in  the  form 
of  an  index  k,  identifying  the  branchpoint  node  in  LIST2(l,k). 
BC(i)  =  0      for  node  i  being  not  a  branchpoint  node 
=  k     for  node  i  (=  LIST2(l,k)  being  a  branch- 
point node. 
6 .   Logical  Operations  On  Matrices 
a.   The  Operation  8 

By  the  S>   -   product  of  two  matrices  A  8  B  in  that 
order  of  the  m  x  p  matrix  A  =  (a- •)  and  the  p  x  n  matrix  B  = 
(bj  .)  is  meant  the  m  x  n  matrix  C  =  (c- ■)  where 

X  J  1  J 
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cij  =  ail  &  blj  •  ai2  &  b2j  e  .  .  .  ffi  air  &  brj 

r 
=  ©   (a-k  &  b,  .),  i=l,2/...m;  j=l/2/...n 
k=l        KJ 

where  &  is  the  logical  operation  AND  with  truth  table: 

a   b    a&b 

0    0     0 

0    10 

10     0 

111 
and  9  is  the  logical  operation  OR  with  the  truth  table: 

a   b    c8b 

0    0     0 

111 

0    11 

10     1 
B.    THE  LOOP  SEARCH 

1 .   Initialization 

Consider  LIST1.   Define  two  temporary  control  vectors 
Tl  and  T2  such  that  for  every  node  there  is  one  binary  position 
reserved.   Then  set  for  every  node  k  present  in  LISTl(l)  the 
element  Tl(k)  to  1  and  for  every  node  m  present  in  LIST1(3) 
the  element  T2 (m)  to  m. 

Tl(k)  =  1   for  LISTl(l,k)  ^  0 

T2  (m)  =  1   for  LISTl(3,m)  ^  0 


20 


Now  the  control  vector  LC  is  initially  set  to 

LC  =  Tl  5  T2 

Example : 

LIST1:    1   a   2 

3  b   2 
2   c   4 

4  d   5 

5  e   6 
graph: 

3 


. 


Fig.  4 

Tl   =  [1   1   1   1   1   0] 

T2   =  [0   1   0   1   1   1] 

LC   =  [1   0   1   0   0   1] 

The  next  step  in  the  initialization  is  the  marking  of 
the  multiple  nodes.   Use  is  made  of  the  control  vector  LIST1(4) 

The  control  vector  LIST1(4)  is  set  to  its  value 
according  to  its  definition.   (See  program  FINDLOOP) 
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2 .   Search  For  Loops 

a.   Finding  Successive  Nodes 

The  Column  1  of  LISTl  is  called  LIST1 (1) ,  the 
Column  2  is  called  LISTl (2).   An  element  in  a'  column  of  LISTl 
will  be  denoted  by  LISTl  (column  number,  element  index) . 

The  search  starts  in  LISTl (1),  using  LISTl (3,1) 
as  the  first  search  key.   The  result  of  this  search  is  an 
index  k  at  which: 

LIST1(3,1)  =  LISTl (l,k) 

Now  the  search  key  is  changed: 
search  key  =  LISTl (3, k) 
and  the  search  continues,  at  the  top  of  LIST(l)  or  at  the 
next  location.   (For  short  lists  it  may  be  reasonable  to  start 
at  the  top  again.) 

When  in  the  continuing  search: 
LISTl(3,k)  =  LISTl(l,i) 
the  search  key  is  again  changed  to: 

search  key  =  LISTl (3, i) 
and  the  search  continues. 

In  a  graph  containing  loops,  LISTl  contains 
multiple  equal  node  numbers  in  both  Column  1  and  2. 

The  question  arises  how  to  proceed  with  the 
search  when  these  multiple  node  numbers  are  encountered.   The 
control  vector  LISTl (4)  contains  information  about  those  nodes 
which  have  leaving  branches  in  more  than  one  direction.   It  is 
necessary  to  follow  all  possible  directed  paths  through  the 
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graph.   LISTl(4)  enables  the  search  routine  to  keep  track  of 
all  nodes  where  to  try  to  follow  also  other  possible  paths. 
If  during  the  search  an  element  LISTl(l,k)  is 
encountered  whose  related  control  element  LISTl(4,k)  =  1,  then 
the  element  LISTl(3,k)  is  selected  as  the  next  search  key  (as 
explained  above),  but  the  control  vector  element  LISTl(4,k) 
is  set  to  "0"  before  the  search  continues. 

b.  Detection  Of  a  Loop 

The  element  LISTl(l,l)  and  all  the  elements 
LISTl(l,n)  which  are  "found"  during  the  search  are  "marked"  in 
a  temporary  control  vector  (the  max.  length  is  equal  to  the 
number  of  nodes) .   Every  time  a  new  search  key  is  selected,  a 
check  is  made  against  the  control  vector.   If  for  the 
specific  node  a  =  LISTl(3,i),  a  "1"  is  found  in  the  control 
vector,  a  loop  is  detected  because  a  closed  directed  path  had 
been  followed,  the  path  being  closed  at  LISTl(3,i). 

c.  Recording  Of  the  Loop 

To  record  the  actual  sequence  of  nodes  in  a  loop, 
which  is  an  orderly  way  of  noting  the  loop  nodes,  it  is 
necessary  to  store  in  a  list  all  successively  "found"  elements 
LISTl(l,n),  starting  with  the  initial  value  LISTl(l,l).   This 
is  done  during  the  search  and  the  list  is  abandoned  if  the 
search  is  terminated  before  a  loop  could  be  found.   But  most 
likely  the  actual  sequence  is  of  no  consequence.   After  the 
jth  loop  detection  the  "loop  matrix"  L  is  updated  by  filling 
the  jth  row  of  L  with  the  appropriate  entries. 
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Also  for  each  node  touched  during  the  search  an 
entry  is  made  in  control  vector  LC . 

d.   Completion  Of  the  Loop  Search 

The  loop  search  is  basically  terminated  if  by 
the  process  described  in  b.  above  no  loop  is  detected.  How- 
ever, before  it  can  be  assumed  that  all  loops  of  the  graph 
have  been  found,  the  control  vector  LC  must  be  checked.   If 
LC  =  1  (contains  only  "1"  's),  all  nodes  in  the  graph  have 
been  checked  for  belonging  to  a  loop.   Otherwise  the  search 
must  be  continued.   For  practical  purposes  the  search  could 
start  at  any  node  which  was  not  yet  included  in  the  search 
and  which  is  identified  by  a  "0"  in  LC. 

The  newly  found  loops  are  recorded  only  if  they 
are  different  from  the  ones  already  found.   This  condition 
can  be  detected  easily,  using  logical  functions  of  the  pro- 
gramming language,  by  comparing  the  new  jtn  entry  row  for 
matrix  L  with  the  other  rows  of  L  already  set. 

The  multiple  starting  of  a  search  is  necessary 
if  the  search  was  started  in  a  part  of  the  graph  which  is 
connected  to  the  rest  of  the  graph  only  by  branches  coming 
from  the  other  part,  as  illustrated  below  in  Fig.  5. 
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Fig.  5 

In  the  graph  given  (Fig.  5)  a  start  at  nodes  5, 
6,  7,  or  8  could  not  lead  to  the  detection  of  the  loops  1-3-4 
and  3-4.   On  the  other  hand  a  start  at  node  2,1,3,  or  4  would 
lead  to  the  detection  of  all  loops  in  one  search. 

Here,  how  fast  the  search  will  be  completed 
depends  on  the  organization  of  the  input  information  by  the 
user. 

C.    THE  PATH  SEARCH 

1.   Initialization 

For  the  purpose  of  finding  all  paths  of  the  graph 
leading  from  a  specified  node  to  another  specified  node,  a 
reorganized  input  list  is  formed.   This  list  is  called  LIST2. 

As  a  first  step  a  multiple  node  control  vector 
LIST2(4)  is  formed,  exactly  analogous  to  LISTl  (4)  . 
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2.   Path  Search 

The  path  search  follows  the  same  logical  sequence  of 
events  as  the  loop  search,  except  that  the  path  search  starts 
not  at  the  element  at  the  top  of  the  list  but  with  LIST2(l,a)/ 
using  LIST2(3,a)  as  the  first  search  key.   The  node  "a"  is 
the  specified  end  node,  to  which  a  path  is  sought. 

The  result  of  the  search  using  the  search  key  is  an 
index  k  at  which 

LIST2(3,a)  =  LIST2(l,k)  . 
Now  the  search  key  is  changed: 

search  key  =  LIST2(3,k) 
and  the  search  continues  as  described  before  (B.2).   This  is 
the  basic  description  of  the  path  search.   The  difference  , 
compared  to  the  loop  search,  is  the  fact  that  here  the  pre- 
decessors of  the  end  node  are  found  and  that  here  all  connec- 
tions are  avoided  which  would  result  in  loops  or  which  would 
not  end  at  the  desired  starting  node. 

a.   Detection  Of  a  Path 

A  path  is  detected  when  it  is  found  that  the 
search  key  is  set  to  an  element  LIST2(3,i)/  where  the  node 
"i"  is  the  specified  starting  node.   The  path  so  detected 
leads  from  node  "a"  to  node  "i",  or,  in  the  notation  used 
here  to  facilitate  computer  programming,  from  LIST2(l,a)  to 
LIST2 (3,i) . 

During  the  search  a  temporary  control  vector  is 
updated  with  "1"  entries  as  described  above  (2.b)  for  the 
loop  detection.   Here  the  element  LIST2(l,a)  and  all  the 
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elements  LIST2(3,n),  which  are  used  as  search  keys  are  marked 
with  a  "1"  in  the  related  place  in  the  control  vector.   If 
for  the  element  LIST2(3,i)  a  "1"  is  already  marked  in  the 
control  vector,  a  loop  is  detected.   This  is  an  unwanted 
situation.   The  loop  must  be  left  by  going  back  to  the  last 
multiple  node.   Then  one  of  the  other  alternate  paths  going 
back  from  that  node  must  be  followed. 

How  is  it  possible  to  keep  track  of  the  multiple 
nodes  passed?   It  must  be  remembered,  that  the  multiple  nodes 
are  marked  in  LIST2(4).   (See  definition  for  the  analogous 
vector  LIST1(4)).   So  every  time  a  new  search  key  is  selected, 
LIST2(4)  will  be  checked,  and  a  "1"  in  the  position  k,  which 
is  the  "found"  -  index,  indicates  that  the  node  LIST2(l,k)  is 
a  multiple  node.   When  passing  these  nodes  during  the  search, 
the  element  LIST2(4,k)  is  set  to  "0"  so  that  it  cannot  be 
found  again. 

The  control  vector  BC  (for  Branch  Control)  is 
used  to  store  all  multiple  nodes  passed  during  the  search, 
probably  best  in  the  form  of  the  index,  where  the  multiple 
node  can  be  found  in  LIST2  (1) .   It  is  then  possible  to 
identify  the  node  to  which  one  must  go  back  to  resume  the 
search,  in  case  a  loop  is  detected,  or  the  path  detected  ends 
at  a  node  not  specified  as  the  end  node. 

The  latter  case  is  possible  and  is  identified 
by  the  fact,  that  using  the  present  search  key,   (the  last  one 
assigned)  no  element  LIST2(l,i)  can  be  found  in  column  LIST2 (1) 
which  is  equal  to  the  search  key. 
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b.   Recording  Of  the  Path 

To  record  the  actual  sequence  of  nodes  in  the 
path,  backwards  from  the  end  to  the  start,  it  is  necessary  to 
store  in  a  temporary  list  all  the  elements  LIST2(3,n)  used  as 
search  keys  during  the  search  and  the  initial  element  LIST2(1, 
a) .   This  list  contains  only  those  nodes  which  actually  form 
the  path.   Intermediate  restarts  of  the  search  are  not 
recorded. 

After  the  search  for  the  jth  path  is  completed, 
the  jtn  row  of  the  matrix  P  is  filled  with  the  appropriate 
entries.   This  can  be  most  practically  accomplished  by 
defining  the  temporary  list  mentioned  above  in  the  form  of  a 
row  of  P.   Then  it  is  only  necessary  to  set  the  j    row  of  P 
equal  to  this  temporary  list  at  the  end  of  the  search. 

D.    EVALUATION  OF  TRANSFER  FUNCTIONS 

To  evaluate  the  transmission  gain  from  an  input  node  "i" 
to  an  output  node  "o"  (here:  a  node,  which  must  have  at  least 
one  incoming  branch)  the  Mason  Gain  formula  is  used: 


G   = 


EAk  Pk 


where  P,  is  the  gain  of  the  kth  direct  path  from  node  "i"  to 

node  "o" 

Ak  is  the  cof actor  of  the  k^  path,  formed  from 

A  by  striking  out  all  terms  containing  loops  which 

are  touched  by  the  kth  path. 
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1 .   Evaluation  Of  the  Graph  Determinant  A  (DELTA) 

If  the  graph  would  contain  only  loops  which  do  not 
touch,  then  the  graph  determinant  would  be  evaluated  as 

DELTA  =  (1-L-l)  (1"L2)  .  .  .  (1-L±) 
where  i  =  nlf  the  number  of  loops  in  the  graph,  and  L-  is  the 
loop  gain,  the  product  of  the  gains  of  the  branches  which  form 
the  loop  L  .  . 

If  touching  loops  are  present,  then  the  products  of 
their  loop  gains  are  eliminated  from  the  expression  for  DELTA, 
a.   Evaluation  Of  the  Matrix  LT 

The  matrix  LT  contains  information  about  which 
loops  touch  each  other. 

The  matrix  LT  is  evaluated  as: 
LT  =  L  8  LT      LT  is  the  transposed  matrix  L, 

S  is  the  operation  defined  in 
I. 6. a 
Proof.   Two  loops  touch  each  other  when  they  have  at  least 
one  node  in  common.   By  definition  of  the  matrix  L  and  the 
operation  8  an  element  LT(k,j)  =  1  if  at  least  one  of  the 
logical  expressions  aik  &  bk .  =  1.   There  are  k  of  these 
expressions,  for  every  node  in  the  graph. 

By  definition  of  the  operation  "  &  "  an  expression 

aik  &  bkj  =  1  onlY  if  aik  =  -1  and  also  bki  =  1'      But  this  is 
the  case  only  if  the  node  k  is  in  loop  i  (for  aik)  and  also  in 

loop  j  (for  bk-).   Consequently  LT(i,j)  =  1  only  if  one  of  the 
graph  is  in  loop  i  and  loop  j  simultaneously.   This  completes 
the  proof. 
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L  = 


It  should  be  noted  that  the  matrix  LT  is 
symmetric  and  has  always  all  elements  on  the  main  diagonal 
equal  to  1 .   As  a  consequence  of  this  all  the  elements  above 
or  below  the  main  diagonal  and  the  elements  on  the  main 
diagonal  could  be  discarded.   The  computer  program,  Mason 
which  is  used  to  evaluate  DELTA  (Program  2)  uses  this  fact. 
For  convenience  the  example  of  Fig.  3  is  given: 
T  1   0   0   0   0  T 
0   0   110   0   0 
0   0   0   1111 
_0   0   0   0   0   1  1_ 
then  LT  is  evaluated: 

T  0 
10  0  0 
0  10  0 
0  110 
0  0  10 
0  0  11 
0   0   11 


LT  = 


110  0  0  0  0 
0  0  110  0  0 
0  0  0  1111 
0   0   0   0   0   11 


1 

0 

0 

0 

0 

1 

1 

0 

0 

1 

1 

1 

0   0   11 


Now  LT  gives  the  result: 
L2  and  L-.  touch, 


L-,  and  L,,  touch 
3       4 


The  number  of  touching  loops  in  the  graph  (of  order  2)  is 
equal  to  the  sum  of  the  elements  of  LT  divided  by  two. 
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b.   Evaluation  of  DELTA  By  a (FORMAC  Program 

A  computer  program,  using  the  language  FORMAC,  is 
used  to  evaluate  the  expression  for  DELTA.  (Program  2) 

The  program  is  only  meant  as  a  demonstration  of 
how  this  problem  of  evaluating  DELTA  can  be  solved  using  an 
algebraic  manipulation  language. 

After  each  multiplication  by  (1-L.)  the  expanded 
result  is  checked  for  the  presence  of  products  of  two  touch- 
ing loops,  and  if  such  a  product  is  found,  it  is  replaced  by 
zero.   It  is  not  necessary  to  check  for  products  of  more  than 
two  touching  loops,  because  larger  products  containing  the 
smaller  one  will  go  to  zero  automatically,  after  the  smaller 
product  is  set  to  zero. 

To  find  the  sets  of  two  touching  loops,  the  upper 
triangular  part  of  LT,  excluding  also  the  main  diagonal,  is 
read  row  by  row.   The  appearance  of  a  "1"  during  that  process 
signals  a  touching  of  the  respective  loops,  indicated  by  the 
row-  and  column  index  of  the  "1"  found.   This  information  is 
stored  in  a  table  in  the  form  of  an  algebraic  expression 
consisting  of  the  product  of  the  two  loops  involved.   From 
this  table  the  expressions  are  taken  when,  after  the  multi- 
plications mentioned  above,  a  check  for  their  presence  in  the 
accumulated  product  is  made. 

2.   Evaluation  Of  the  kth  Cof actor  AR  (DELTK) 

Having  evaluated  DELTA,  the  graph  determinant,  the 
ktn  cof actor  is  evaluated  in  a  similar  way,  using  the  same 
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computer  program.   Instead  of  the  matrix  LT  now  the  path- 
loop  touching  matrix  PL  is  used  as  input. 

The  matrix  PL  is  evaluated  as: 

PL  =  L  S  PT 
where  PT  is  the  transposed  path  matrix  P.   The  resultant  mat- 
rix PL  shows  the  path  -  loop  connections.   If  the  matrix  is 
read  row-wise,  then  for  each  row  i  and  column  j  it  is  indicated 
by  a  "1",  if  loop  i  touches  path  j.   The  proof  follows  the 
proof  for  the  case  of  the  touching  loop  matrix  LT. 

For  the  example  given  for  the  loop  case:   (Fig.  3) 

Til 

10    0 
PL  = 

111 

111 
which  indicates,  that  only  loop  2  is  not  touched  by  path  2 
or  path  3 . 

Now  two  possible  ways  of  solution  can  be  followed. 
One  way  is  to  analyze  the  matrix  PL  in  a  way  similar  to  the 
one  used  for  the  matrix  LT  before.   For  each  column  of  PL  a 
specific  cof actor  DELTK,  is  evaluated.   For  a  "1"  in  row  i, 
loop  L(i)  is  eliminated  from  the  previously  evaluated  expres- 
sion for  DELTA.   Now  the  simple  replacement  of  L(i)  by  zero 
results  in  elimination  of  all  terms  in  DELTA  which  contain 
that  loop  L(i).   After  processing  in  that  way  all  columns  of 
PL,  all  cof actors  are  evaluated.   Now  DELTA  and  the  evaluated 
cofactors  DELTKk  for  the  k  paths  present  can  be  used  in  the 
final  evaluation. 
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The  other  possible  solution  is  first  to  use  DELTA  in 
the  final  evaluation  by  substituting  for  every  Li  in  DELTA 
the  actual  algebraic  expression  for  L±   and  then  evaluate  the 
cofactors  new  from  the  original  expression  involving  the 
accumulated  products  of  (1-L. ) ,  similar  to  the  evaluation  of 
DELTA . 

3.   Evaluation  Of  the  kth  Path  Gain  (PK) 

The  algebraic  evaluation  of  the  path  gain  is  best 
done  immediately  after  a  path  is  found  (see  II.C.2.b).   One 
method  of  solution  would  be  to  store  during  the  search  for 
the  path  also  the  index  k  of  the  different  search  keys.   This 
index  sequence  could  then  be  used  in  the  algebraic  part  of 
the  program  to  find  those  branch  gains  in  the  input  table 
LIST1  which  must  be  multiplied  to  form  the  path  gain. 
Two  basic  approaches  are  possible: 

-  the  symbols  stored  in  LIST1(2)  represent 
larger  and  more  complex  algebraic  expressions. 
Then  they  constitute  names  or,  to  use  FORMAC 
vocabulary,  they  are  FORMAC  variables.   In  this 
case  there  must  be  an  assignment  of  the  alge- 
braic expression  to  that  variable  somewhere  in 
the  program.   (See  example  below) 

-  the  symbols  stored  in  LIST1(2)   are  the  actual 
algebraic  expressions  which  constitute  the 
branch  gains.   In  FORMAC  vocabulary  they  are 
called  FORMAC  expressions.   In  this  case  there 
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Fig.    6 


m 
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must  be  an  evaluation  which  calculated  the  path 
gain  directly.   (See  example  below) 

Example. 

Case  1. 

LISTl(l)        LIST1(2)         LISTK3) 


1  b  2 

2  a  1 
2  c  3 
2  k  6 

2  1  7 

3  e  4 

4  d  3 

4  g  5 

5  h  6 

6  i  7 

7  f  4 
7               j  6 

Note:   the  elements  in  LIST1(3)  must  be  character  type. 

In  PL/1  this  is  indicated  by  writing  for  the  character 
b  the  expression  'b1.   (See  PL/1  reference  manual) 

A  FORMAC  statement  like  the  following  would  produce 
the  path  gain  in  terms  of  the  letters  in  LIST1(3): 

LET(PK(1)  =  "LIST1(2,1)"  *  "LIST1 (2 , 3) " 

*  *  LISTl(2,i)  ) ; 

The  following  statement  would  produce  this  output: 

PRINT_Out (PK(1) ) ; 

PK(1)  =  B  C  E  ...  J 
Because  generally  it  would  not  be  known  explicitly  which 
symbols  out  of  LIST1(2)  form  the  gain  PK,  the  following 
assign  procedure  would  be  used: 

At  the  time  a  more  specific  evaluation  in  terms  of 
the  actual  algebraic  expressions  has  to  be  produced,  an  input 
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statement  is  executed.   This  statement  would  cause  the  read-in 
of  the  actual  algebraic  expressions.   Then  a  replacement 
would  take  place,  in  which  the  letters  a,b,...   are  replaced 
in  the  gain  expression  PK(1)  =  B  C  E  G  ...  F  by  the  algebraic 
expressions.   Then  a  further  expansion  of  the  gain  expression 
could  take  place. 

On   data  cards : 

iK7.  ,3  *  A  *  *  2-  'i/s«  

Then  the  statements: 

DECLARE  EXPR  CHARACTER ( 80 )  VARYING; 

to  define  EXPR  a  character  variable, 
to  define  an  index. 

K  =  0; 
A:    GET  LIST(EXPR);  K  =  K+l; 

ON  ENDFILE(SYSIN)  GO  TO  END; 

LET(PK(1)  =  REPLACE (PK (1) ,  "LISTl (2 ,K) " ,  "EXPR")  ); 

GO  TO  A; 
END:  LET(PK(1)  =  EXPAND  (PK(1))  ); 

PRINT_OUT(PK(l) ) ; 
The  printout  would  read: 

PK(1)  =  K7  3  A2  /  S 
In  this  case  the  further  expansion  has  no  effect, 
because  there  are  no  expansions  to  be  performed. 
Now  consider  case  2 : 

LISTl  (1)  LISTl (2)  LISTl (3) 

as  K7  as 

before  3*A**2  before 

1/S 
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In  this  case  there  are  no  data  cards.   The  statements: 

Note:  NB  =  number  o-  branches 

in  a  path 

LET(PK(1)  =  1) ; 

DO  K  =  1  TO  NB; 

LET(PK(1)  =  PK.(1)  *  "LIST1  (2,INX(K)  )  "  ); 

END; 

PRINT_0UT(PK(1)  )  ; 

Note:   INX  is  an  array  of  indices 

which  indicate  the  branch 

to  be  included  in  the 

product. 

In  the  example  given  the  array  INX  would  contain  for 

the  first  path  PK(1)  the  following  indices: 

INX(l)  =  1 

INX (2)  =  3 

INX  (3)  =  6 

INX  (4)  =  8 

INX (5)  =  9 

INX(6)  =  10 

and  NB  =  6 . 

A  similar  evaluation  has  to  take  place  for  the  loop 
gains,  the  values  of  DELTA  and  DELTK. 

4 .   Final  Evaluation  Of  the  Transmission  Gain 

The  transmission  gain  G  can  finally  be  evaluated 
using  a  FORMAC  program.   Depending  on  the  complexity  of  the 
problem  first  G  could  be  evaluated  in  terms  of  the  branch  gains, 
where  the  branch  gains  are  represented  by  their  variable  names, 
and  by  the  loop  gains.   In  this  case  G  would  be  an  algebraic 
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expression  such  as 

G  =   (L(i)L(k)  ..  .)  *  (A,B,C.) 
(L(a)L(b)  ...) 

That  is,  DELTA  and  DELTK(i)  would  be  in  terms  of  the 
loop  gains;  L(i)  and  PK(i)  would  be  in  terms  of  the  branch 
gains  a,b,c. . . 

Then  after  G  is  evaluated  as  above,  a  substitution 
could  be  made,  where,  again  depending  on  the  complexity 
involved  in  possible  applications,  several  methods  are  possible, 
In  this  process  the  FORMAC  functions  REPLACE,  CHAIN,  EVAL  and 
EXPAND[1]  would  probably  have  to  be  used.   Generally  a  com- 
promise has  to  be  made  between  using  more  time  for  execution 
and  less  core  space  and  vice  versa. 
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III.   DIRECT  EVALUATION  OF 
DETERMINANTS  USING  FORMAC 

In  this  chapter  the  direct  evaluation  of  determinants 
using  FORMAC  programs  is  shown.   By  "direct"  evaluation  is 
meant  a  process  of  evaluation  where  the  properties  of  deter- 
minants are  used  rather  than  graph  theory  as  used  in  chapter 
II  for  the  evaluation  of  DELTA  (which  is  a  determinant) . 

A.    A  FORMAC  PROGRAM  USING  THE  BASIC  DEFINITION  OF  THE 
DETERMINANT 

In  numerical  mathematics  a  determinant  is  normally  eval- 
uated by  normalizing  the  rows  or  columns  and  using  elementary 
row-  or  column  operations  to  eliminate  all  elements  off  the 
main  diagonal.   Then  the  value  of  the  determinant  is  the 
product  of  the  elements  on  the  main  diagonal.   If  this  process 
would  be  used  in  the  evaluation  of  determinants  with  algebraic 
terms  as  elements,  the  generation  of  rational  functions  or 
otherwise  more  complex  algebraic  expressions  would  create 
additional  complexity.   The  handling  of  rational  polynomial 
expressions,  when  the  polynomials  are  in  symbolic  form,  is  very 
time  consuming. 

The  program  DTERM  uses  the  basic  definition  of  the 
determinant : 

If  A  is  a  square  matrix  of  order  n,  then  the  determinant 

Of  A  is  DET(A) : 

t 
DET(A)  =  E(-l)  a±      a2  -  "an 

3         jl   32      jn 
where  the  second  subscripts  j  assume  all  possible  arrangements 

in  which  each  column  is  represented  exactly  once  in  each  term 
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of  the  sum,  and  the  exponent,  t,  is  the  number  of  inter- 
changes necessary  to  bring  the  second  subscripts  into  natural 
order  (that  is,  1,2,3,  ...,  n). 

In  other  words,  all  possible  permutations  of  the  numbers 
1,2,3,  ...,  n  form  the  complete  set  of  second  subscripts 
mentioned  above. 

If  the  number  t  is  even,  the  permutation  is  called  "even", 
otherwise  "odd". 

The  permutation  of  a  set  of  numbers  soon  exceeds  the 
time  available  for  execution  of  any  program.   However,  for 
n  =  2  to  n  =  9  one  may  find  it  possible  to  use  this  method. 
The  algorithm  used  in  DTERM  to  find  all  the  permutations  of 
the  numbers  1,2,3,  ...,  n,  for  different  n  has  the  following 


execution  time: 

n  =  4 
n  =  5 
n  =  6 
n  =  7 
n  =  8 
n  =  9 


0.05  sec 
0.09  sec 
0.31  sec 
1.95  sec 
15.12  sec 
137.95  sec 

In  the  example  program  DTERM  the  permutation  of  column 
indices  is  only  used  to  demonstrate  the  algebraic  part  of  the 
program.   Normally  this  algebraic  part  would  be  used  with 
another  method  for  finding  the  elements  of  the  matrix  A  which 
have  to  be  multiplied  to  form  a  term  in  the  expansion  of  the 
determinant.   For   very  large  matrices,  (for  example,  a  37  x 
37  matrix)  the  graph  analysis  as  described  before  would 
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identify  those  elements.   Also  the  method  used  by  Lipson  [4] 

could  possibly  be  used  to  identify  those  elements.   How  much 

computational  effort  each  method  requires  in  terms  of  time  and 

space  has  to  be  determined. 

1.   General  Description  Of  the  Evaluation  Of  the 
Determinant  By  the  FORMAC  Program  DETERM 

a.    The  Non-Algebraic  Part  Of  DETERM 

The  first  step  after  the  initialization  is  the 
evaluation  of  a  permutation  of  the  numbers  1,2  ...,  n.   The 
program  used  for  that  purpose  is  the  subroutine  PERM1 .   The 
program  PERM  produces  upon  each  call  another  distinct  permuta- 
tion of  the  n  !   possible  permutations  in  such  a  way  that  the 
permutations  it  generates  are  alternately  odd  and  even,  a 
property  the  author  of  algorithm  115  PERM  described  as  "prob- 
ably useless".   However,  this  property  enables  us  to  determine 
the  sign  of  the  permutation. 

Having  found  a  specific  permutation  it  is  used  as 
the  set  of  indices  jl,  j2,  ...,  jn  mentioned  in  the  definition 
of  the  determinant.   At  the  same  time  the  sign  of  the  product 

a,   a^   ...a    is  known. 
1.   2.      n- 

Jl   J2      Jn 

The  elements  a.   in  the  product  are  then  tested 
1j 

for  being  zero.   If  a  zero  element  is  found,  a  new  permutation 
is  produced  and  the  process  starts  again,  because  further 
algebraic  manipulation  would  only  produce  a  zero  term. 


xThe  authors  conversion  of  the  algorithm  115  PERM,  H.  F. 
Trotter,  Comm.  ACM,  Aug.  1962  from  ALGOL  into  PL1. 
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In  case  the  elements  are  all  non-zero  the  product 
is  algebraically  evaluated.   This  is  a  complicated  step, 
which  will  be  described  in  the  next  section. 

Then  the  next  non-zero  term  (a  term  is  one  of 
the  products  aijiakjk- • -an jn)  is  evaluated  in  the  same  manner. 

Finally  all  coefficients  of  equal  powers  of  the 
polynomial  variable  in  the  different  terms  of  the  determinant 
are  collected. 

b.   The  Algebraic  Part  of  DETERM 

Most  of  the  algebraic  manipulation  statements 

are  detectable  by  the  word  "LET  ( )  ;  and  the 

embracing  brackets.   The  semicolon  is  a  feature  of  PL/1;  it 
signals  the  end  of  a  statement. 

Consider  now  two  polynomials  in  the  variable  S: 
PI  =  (A  +  B2)2  s  +  (C  +  D(E  +  F)3  S2 
P2  =  (A  +  C)3   S  +  B  S3  +  F(G-H)   S5 
If  the  product  Pl*P2  is  formed  and  all  the  terms  are  completely 
expanded,  a  new  polynomial  is  formed  which  has  90  terms  before 
simplification. 

If  the  same  product  is  formed  without  expanding 
the  coefficients  of  the  variable  S  the  resultant  polynomial 
has  only  6  terms.   Each  of  the  coefficients  of  the  variable  S 
is  in  unexpanded  form. 
Example. 

PI  =  (a  +  b)  s  +  (c  +  d)  s2 
P2  =  (e  +  f)  s  +  (g  +  h) 
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The  expansion  used  in  DTERM: 

PI  P2  =  (a  +  b)  (e  +  f )  s2  +  (a  +  b)  (g  +  h)  s 

+  (c  +  d)  (e  +  f )  s3  +  (c  +  d)  (g  +  h)  s2 
Then  this  result  is  simplified  to 

PI  P2  =  (a  +  b) (g  +  h)  s  + 

(  (a  +  b)  (e  +  f)  +  (c  +  d)  (g  +  h)  )  s2  + 
(c  +  d)  (e  +  f)  s3 
or  =  A  s  +  B  s2  +  C  s3 

where  the  coefficients  A,  B,  C  ...  are  stored  in  this  form  on 
disks.   After  all  terms  of  the  determinant  are  expanded  in 
this  manner,  coefficients  of  equal  powers  of  s  are  read  back 
into  main  core  and  are  expanded  fully  and  added.   They  are 
then  printed  and  the  space  used  by  them  is  subsequently  made 
available  to  other  uses. 

The  latter  method  is  used  in  the  first  part  of 
the  algebraic  manipulation  to  save  space.   A  main  limitation 
of  FORMAC  (probably  of  any  algebraic  manipulation  language)  is 
the  enormous  space  requirement  for  certain  expansions;  there- 
fore to  make  algebraic  manipulation  of  large  expressions 
possible  at  all,  it  must  be  done  at  the  expense  of  time.   It. 
has  been  shown  [6]  that  running  FORMAC  programs  can  be  made 
much  faster  after  some  refinement;  often  the  omission   of  the 
use  of  one  FORMAC  function  showed  a  considerable  increase  in 
execution  time.   The  reason  is  that  each  FORMAC  function  use 
constitutes  really  the  use  of  a  subroutine. 
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The  main  emphasis  in  DETERM  was  placed  on  the 
possibility  of  evaluating  very  large  determinants,  the 
refinement  was  secondary. 

After  all  n  factors,  where  each  is  considered 
being  a  polynomial,  have  been  multiplied,  the  resultant 
product,  which  is  called  "PD"  in  the  program,  is  further 
processed. 

Each  coefficient  in  PD  of  the  variable  S  and  the 
constant  coefficient  is  identified  and  stored  separately  on 
disk  file.   Then  the  evaluation  of  a  new  term  PD  starts. 

After  all  terms  PD,  which  are  the  products 
mentioned  in  the  definition  of  the  matrix,  have  been  processed 
in  the  manner  described  above,  the  coefficients  stored  on  disk 
are  now  retrieved  and  added  up,  always  adding  coefficients  of 
the  same  power  of  the  variable  S  originating  from  the  different 
terms  of  the  determinant.   The  coefficients  are  designated  as 
follows : 

COEF(i)  =  C(l,i)  +  C(2,i)  +  ...  +  C(k,i) 
where  i  is  the  power  of  S  which  is  under  consideration  and  the 
numbers  1,2,  . . .k  are  the  numbers  of  the  k  non-zero  terms  in 
the  determinant. 

The  coefficients  C,  •  are  the  ones  which  are 
stored  on  disk.   The  final  coefficient  COEFi  is  the  coefficient 
of  the  variable  power  S1  in  the  expanded  determinant. 
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Example  for  using  DTERM2 

System  equations: 

e.  =  k  .  0. 
i    pi   1 

e   =  k    0 

0  po  o 

e   =  e .  -  e 
i    o 

e,.  =  k   e 
f    a   a 

e^  =  h  0 
t      t 

e  =  e   -  e 
act 

Xm  =  kt  *f  \ 

Ec(s>    c  s  +  a 

E  (s)  ■  uc    s  +  b 

The  state  equations  are  formed  with  the  following  variables  as 

state  variables: 

xn  =  e 

1  o 

X3  "  Lf 

X  =  dummy  variable 

Ec(s)    X4   Ec     1       s  +  a 
Gc(s)  =  E(s)   =  E~  XT    s  +  b 


2Example  taken  from  notes  course  3411  Naval  Postgraduate 
School. 
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which  gives: 


s  X,  =  -  (b)  X,  +  e 
4  4 


e   =  X,  +  a  X, 
c     4      4 


The  equations  in  the  S  -  domain: 

(1)  s  xx  -  x2  =0 

(2)  (jo(p1  P2)2  +  Jm)  s  X2  +  (fm(p1  P2)2  +  fm)  X2 

"  (pl  P2  kt  ia)  X3  =  °       J 

(3)  (ka  kpo}  Xl  +  (ka  h  VP2)  X2  +  (rf>  X3 

+  (L  s)  X,  +  (k  (b  -  a) )  X.  =  k  k  .0. 

3    v  a  4  apii 

(4)  (s)  X,  +  (k   )  Xn  +  (a  +  b)  X,  =  k  .0. 

4     po    1  4  pi  i 

The  equations  in  matrix  form: 

A  X  =  B 
where  A  is  a  matrix  of  polynomial  coefficients 

X  is  a  matrix  of  variables 

B  is  a  matrix  of  constants 
The  elements  of  A  as  used  as  input  to  DTERM: 

An  =  s 

A12  -  -  1 

2  2 

A22  =  (J°(Px  p2^   +  Jm)  s  +  (fm(Pi  p2^   +  fm^ 

A23  "  "  (Pi  ^2  kt  ia) 
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A31  =  (ka  V 
A32  =  (ka  h  l/p2) 

A33  =  (L  s  +  rf) 

A...  =  (k   ) 

34      po 

A44  =  (s  +  a  +  b) 

The  program  DTERM  with  the  matrix  A  as  input  and  the  output 
is  contained  in  section  COMPUTER  PROGRAMS. 

Figure  8  shows  the  physical  system,  Fig.  9  shows 
the  Mason  flow  graph. 
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Figure    8 
48 


Figure    9 
49 


2.   Discussion  Of  the  Flow  Graph  For  DETERM 

The  flow  graph  for  DETERM  is  discussed  here,  to  some 
limited  extent. 

For  clarification  the  definition  of  the  determinant 
is  now  again  given  in  terms  of  the  symbols  used  in  the  program. 

DET(A)  =  PD1  +  PD2+  .  .  .  +  PDin 
where 

PD.  =  C(j,i)  S1       i  =  0,1,2 

and  finally 

DET(A)  =  COEF(0)  +  COEF(l)  S  +  COEF(2)  S2  +  ... 

,MXPOW 


+  COEF(MXPOW)  S1 


where 


COEF(i)  =  C(l,i)  +  C(2,i)  +  ...  +  C(IN,i) 
a.   Data  Design  For  DTERM 

Consider  now  the  computer  program  printout. 
(Program  1) .  All  statements  beginning  with  "DCL"  are  declara- 
tions, they  define  those  data  elements  used  in  the  program 
which  are  defined  explicitly.   Some  of  these  definitions  will 
now  be  discussed. 

IP  is  an  integer  array.   It  contains  after  each 
call  of  PERM  a  distinct  permutation  of  the  numbers  l,2,...n. 
These  numbers  are  used  as  column  indices  in  designating  elements 
of  the  input  matrix  A.  The  array  A  is  a  matrix  of  characters. 
Each  element  can  have  a  varying  length  of  maximal  75  characters. 

PDD  is  a  character  type  variable.   It  contains  at 
certain  moments  an  element  of  the  matrix  A.Nl  is  the  order  of 
the  input  matrix  A. 

50 


PRRT  3 


YES 


Pfiftr  6 


General   flowgroph  of   DTERM 
(   DTERM    J 


Initialization 


I 


PRRT  I 


PPiRT  1 


PERM 


find  a  permutation  J 
of  column  indicts 
PRRT  B. 


Evaluate  the  product  PD. 
The  coefficients  C;  are  not 


expanded. 


Extract  the  coefficients  C 


from  PD:.  Store  on  disk  us- 
ing SAVE. 


PBRT  ¥ 


PRRT  S 


Add  all  coefficients  of  equal 
powers  of  the  variable  s. 
This  results  in  C  (i). 


PR*T  7 


Print  result 
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SPECIALIZED  FLOWCHART   OF  DTERM 
(    DTERM    \ 


Increase 
Index   I 


Index  J 


I 


Definition  of  data  elements 


Initialization  and  input  of  matrix  A 


5 


PRKT  1 


PH*T  Z 


ndex 


Determination  of  sign  of  the 
permutation  i 

i 
PERM 

Determination  of  permutation  i 


YES 


I 


PR*r  3 


ncrease  index  IN 


IN  Counts  ttto  number  of 
nonioro  products  in  tnv 
determinant 


Initialize  for  product  accu 

mulation  of  PD,^ 
m 


A 


1. 


Index  M 


I 


Use  the  jth  argument  of  new 
element  as  multiplier. 


5 


PRRT  r 


Increase  index 


I 


Multiply  the  Ith  argument 
of  the  previous  result  P5 
by  multiplier 


i 
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€> 


PDHTf, 


1 


Result  of  this  step  is  assigned 
another  name  to  make  it  the 
'next  'previous*  result.  PD*P5 


NO 


Find  highest  power  of  S  in  all 

products  PD      ■  MAXPOW. 
IN 


Find  all  coefficients  C  ( i )  of  s' 
PD      .  Store  on  disk  using  SAVE. 


PRRT  f 


PRRT   S 


NO 


■© 


Index  M 


Take  coefficient  C     (i)  of  s'in  PD 

M  IN 


6 


d> 


NO 


Pfitr  7 


M  is  a  running  index 
from  7  to  IN 
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1 


Dtttrmint  the  exponent  of  the 
variable  s  for  that  result. 


PRRT  ¥ 


res 


Find  coeffient  of 
same  power  in  pre- 
vious result  P5. 


res 


Add  result  of 
multiplication 
to  previous 
result  P5. 


Store  P5  on 
disk  using 
SAVE. 


Add  result  of  multiplication  to 
previous  result   P5  . 


Build  an  argument 
list  for  'REPLACE' 
function. 


Store  P5  on  disk 
using  SAVE. 


Replace  the  old  coefficient  of 

pow 
s  with  the  sum  of  the  old 

one  +  result  of  multiplication. 


T>RRT  ¥ 


G> 


NO 


YES 


COEF(i)  =  CM(i) 


COEFli)    =  prtvious  sum 
+  ntw  C     (i) 


NO 


YES 


PR*r  7 
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b.   The  Logic  Sequence   (See  flow  chart,  Fig.  9) 
The  initialization  includes  some  statements 
which  are  self  explanatory.   Some  will  be  discussed  now.   N2 
is  used  in  a  loop  to  fill  the  array  IP  initially  with  the 
sequence  of  numbers  1,2,  fffn  in  natural  order,  an  input 
requirement  for  the  subroutine  PERM.   Also  the  binary  variable 
"FIRST"  is  used  for  control  of  PERM.   The  following  loop, 
using  the  index  1=1  to  N!  is  the  main  loop  of  the  program.   In 
that  loop  all  n  !   possible  products  of  the  determinant  are 
investigated.   Then  the  sign  of  the  permutation  is  determined. 

The  loop  M=l  to  N=l  includes  the  zero  check  for 
the  itn  product  of  n  !   products. 

The  next  loop  M=l  to  Nl  is  reached  only  if  all 
factors  in  a  product  are  non-zero.   Then  the  algebraic  manipu- 
lation begins.   The  procedure  was  described  in  2.     before 
The  action  of  the  various  FORMAC  functions  (NARGS,  ARG,  LOWPOW, 
ATOMIZE,  etc.)  is  described  in  the  FORMAC  reference  manual  [6]  . 

B.    EVALUATION  OF  A  DETERMINANT  USING  BARREISS  FRACTION-FREE 
ALGORITHM 

Another  method  of  evaluating  a  determinant  is  using  an 
algorithm  given  by  Lipson.   This  algorithm  is  based  on  a  two- 
step  variant  by  Barreiss  of  a  fraction-free  elimination  algo- 
rithm, which  in  turn  is  a  modification,  attributed  to  Jordan, 
of  Gaussian  elimination  [4]  . 

Lipson  uses  the  algorithm  to  evaluate  a  system  of  linear 
algebraic  equations  by  transforming  a  matrix  into  triangular 
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form  and  following  back-substitution.   The  whole  process  is 
fraction-free . 

The  linear  equations  processed  have  symbolic  coefficients 
in  the  form  of  multivariate  polynomials. 

Lipson  has  also  programmed  the  algorithm  in  the  language 
"SCOPE  FORMAX"  [4] . 

1 .   The  Fraction-Free  Gaussian  Elimination 

Consider  the  n  x  (n+m)  augmented  matrix  of  a  system 
of  linear  equations  Ax  =  B : 

Aij(0)  =  (aL.   (0))  =  (A|B) 

Then  the  modified  algorithm  is  given 

_  (k-1)  a(k-l)  _  a(k-l)   -,(k-l) 

(k)     \k    aij       akj      aik 

aij        a(k_2) 

k-1, k-1 

where 

k  =  1,2,...,  n-1 

i  =  k+1 , . . . ,  n 

j  =  k+1 , . . . ,  n+m 

and  with  A^  '  as  given  before,  and  where 

(-1) 
a00    ~  X 

=  (k-l)         .      .      .    . 

a .         i  =  1 , . . . , k ;    1=1,..., n+m 

a(k)   _    ^ 

0  i  =  k+l,...,n;  j  =  l,...,k 

Also,  given  here  without  proof  (for  the  proof  see 
Ref.  4)  : 
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The  determinant  of  the  matrix  A 

DET(A)  =  a(n_1) 
nn 

The  formula  given  above  was  programmed  using 
normal  FORMAC,  also  ignoring  for  purpose  of  clarification  any 
pivot  element  check  for  zero. 

The  algorithm  alone  does  not  give  a  fraction-free 
algebraic  representation.   Some  processing  is  necessary  to 
bring  dividend  and  divisor  into  a  form  such  that  cancellation 
takes  place. 

A  very  useful  method  may  be  to  use  this  algorithm 
to  compute  the  determinant  of  a  matrix  with  single  letters  as 
elements,  followed  by  a  substitution  process,  where  more 
complicated  expressions  are  substituted  for  the  letters. 
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IV.   APPLICATION  OF  FORMAC 
IN  OTHER  ALGEBRAIC  PROBLEMS  OF  CONTROL  THEORY 

In  this  chapter  some  other  applications  of  the  algebraic 
manipulation  language  FORMAC  in  control  theory  are  proposed. 
Some  examples  will  be  given,  the  programs  being  more  of  demon- 
stration programs  than  operational  programs. 

The  reason  for  this  is  mainly  the  lack  of  time  which  is 
needed  to  prepare  an  optimal  program,  optimal  with  respect  to 
execution  time  and  space  requirement. 

A.    ROOT  LOCUS  THEORY 

The  impulse  response  of  a  linear  time-invariant  system  is 
determined  by  a  sum  of  exponential  time  functions.   The  expo- 
nents of  these  terms  are  the  roots  of  the  characteristic 
equation  of  the  system. 

A  system  is  stable  if  and  only  if  none  of  the  n  roots  of 
the  characteristic  equation  has  positive  real  parts. 

Consider  the  polynomial  in  the  Laplace-variable  S  which 
represents  a  characteristic  equation: 

n      k 
F(s)  =   I    aksK  =  0 

k=o 
The  coefficients  a,  can  be  complex.   Some  of  the  coeffi- 
cients can  include  a  design  parameter,  K. 

The  polynomial  F(s)  can  be  partitioned  into  two  parts  as 

F(s)  +  K  B(s)  (1) 
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The  complex  variable  S  is  replaced  as  s  =  a  +  j  go.   We 
require  the  real  and  imaginary  parts  of  F(  a  ,  u  )  to  be  zero 
independently : 

Re  A(a,oj)  +  K  Re  B(a,co)  =  0  (2a) 

Im  A(a,u>)  +  K  Im  B(a,oo)  =  0  (2b) 

Using  the  set  of  simultaneous  equations  (2a)  and  (2b)  the 
parameter  K  can  be  eliminated: 

Re  B  *  Im  A  -  Re  A  *  Im  B  =  0  (3) 

This  equation  is  the  equation  of  the  root  locus;  that  is, 
every  root  of  the  characteristic  equation  F(s)  =  0  must  satisfy 
the  equation  (3) . 

In  the  computer  program  (Program  3)  this  equation  is 
called  RLOC  =  f(a,u)).   RLOC  contains  in  every  term  the 
variable  co,  because  the  real  axis  is  always  a  locus  for  the 
roots.   This  factor  of  w  can  be  removed  from  RLOC. 

If  now  a  value  for  one  of  the  variables  a  or  co  is 
assumed  RLOC  will  give  those  points  in  the  s-plane  which 
satisfy  equation  (3).   These  points  are  on  the  root  locus. 
1.   A  FORMAC  Program  For  the  Root  Locus 

As  an  example  for  the  application  of  FORMAC  a  program 
is  given  (Program  3). 

The  program  reads  in  a  polynomial  CHEQ  and  the  para- 
meter PARAM  in  algebraic  form.   Then  the  polynomial  is  expand- 
ed and  then,  following  the  development  as  described  before, 
the  equation  (3)  of  the  root  locus  is  derived  in  algebraic 
form. 
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At  this  point  values  for  to  or  a  in  the  range  of 
interest  can  be  chosen  and  the  equation  for  the  root  locus 
can  be  evaluated,  giving  values  for  10  or  o  respectively. 

Normally  the  expression  in  the  remaining  unknown 
will  not  be  linear.   In  that  case  the  expression  is  treated 
as  a  polynomial,  the  roots  of  which  can  be  found  using  a 
root-finding  routine. 

Finally  the  results  can  be  plotted,  giving  the 
complete  root  locus  in  the  range  of  interest. 

To  evaluate  the  values  of  the  parameter  at  points 
along  the  root  locus  it  is  necessary  to  substitute  back  the 
values  of  a  and  co  for  those  points  into  equation  (2)  . 
2 .   Discussion 

The  conventional  root-locus  procedure  [9]  is  basically 
a  trial-and-error  method: 

For  a  range  of  parameter  values  K^,K2/.../  K  ,  for 
which  one  suspects  an  acceptable  area  of  root  locations  in 
the  s-plane,  the  roots  of  the  characteristic  equation  are 
found  by  a  root-finding  method,  probably  a  computer-programmed 
routine,  using  a  standard  method. 

The  the  root  locus  is  inspected  in  the  area  of 
interest  for  the  values  of  a  and  co  and  the  corresponding  para- 
meter values  are  selected. 

The  problem  is  that  values  of  a  quantity  which  are 
sought  have  to  be  assumed  beforehand. 

The  method  proposed  here  [10]  avoids  this  sequence. 
The  equation  of  the  root  locus  is  evaluated  in  the  range  of 
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interest  for  the  variables  a  and  w.   Then  by  back-substitution 
of  a  and  oa  into  the  given  characteristic  equation  the  para- 
meter values  are  calculated.   Another  advantage  is  the  fact 
that  the  degree  of  the  root  locus  equation  (3)  is  generally  of 
lesser  degree  than  the  characteristic  equation.   This  face 
results  from  possible  cancellations  of  terms  of  equal  powers 
in  the  process  of  parameter  elimination.   In  example  1  the 
degree  of  the  characteristic  equation  is  4  while  the  equation 
of  the  root  locus  is  of  degree  3  in  a  and  of  degree  2  in  w . 

B.  FREQUENCY  FUNCTIONS 

Consider  a  transfer  function  P(s)  =  Y(s)+X(s).   P(s) 
might  contain  some  parameter  in  symbolic  form.   It  is  possible 
to  use  a  FORMAC  program  to  evaluate  the  absolute  value  of  P(s) 
and  the  angle  of  P(s)  as  expressions  containing  the  parameter 
in  symbolic  form.   For  complicated  transfer  functions  it 
might  be  easier  to  evaluate  the  complex  algebra  only  once  and 
then  vary  the  parameter  to  obtain  values  for  magnitude  and 
angle  of  the  transfer  function. 

For  the  latter  purpose  it  is  mentioned  that  in  FORMAC 
transcendental  functions  with  symbolic  arguments  can  be  used. 

C.  LAPLACE  TRANSFORMS  AND  INVERSION  OF  LAPLACE  TRANSFORMS 
Consider  a  differential  equation  in  the  time  domain.   If 

such  an  equation  is  given  as  input  to  a  FORMAC  program  it  is 
possible  to  perform  a  Laplace  transform  on  the  equation.   The 
basic  problem  to  be  solved  would  be  the  replacement  of  one 
operator  by  another  one.   If  both  sides  of  the  differential 
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equation  would  be  read  in  as  two  separate  FORMAC  variables,  a 
transfer  function  could  be  formed. 

Another  definite  possibility  is  the  analysis  of  an  alge- 
braic expression  in  the  variable  s  and  subsequent  recognition 
of  certain  algebraic  patterns,  which  can  then  be  replaced  by 
an  expression  in  the  time  domain.   This  method  would  be 
similar  to  the  "manual"  approach,  because  one  often  solves 
these  problems  in  the  same  way:  recognizing  the  transform 
pair  from  experience  (or  tables) . 
D.    SENSITIVITY  ANALYSIS 

A  very  powerful  feature  of  FORMAC  is  its  ability  to  form 
partial  derivatives  of  user  specified  order.   Also  the  mathe- 
matical functions  such  as  sine,  cosine,  exponential,  logarithm 
can  be  used  with  symbolic  arguments. 

Consider  a  mathematical  model  T(k),  which  might  be  a 
transfer  function  or  frequency  response  function,  of  a  linear 
time-invariant  system. 

The  sensitivity  of  T(k),  k  being  a  variable,  with  respect 
to  k  is  defined:  [9] 

sT(k)      d  In  T(k)   =   d  T(k)    k 


k         d  In  k        d  k    T(kT 
The  sensitivity  of  the  phase  angle  of  T(k),  $T  with 
respect  to  parameter  k  is  defined: 

$T       d  <t>T  k 
Sk         d  k  $T 

A  FORMAC  program  which  would  find  the  derivatives  involved  in 
the  sensitivity  analysis  could  be  time  saving,  if  symbolic 
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arguments  were  involved  and  one  would  like  to  vary  the  para- 
meter.  A  program  for  the  finding  of  the  derivatives  would 
only  be  used  once. 
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V.   CONCLUSION 

A  method  of  solution  of  a  pattern  recognition  problem  in 
the  analysis  of  a  Mason  flow  graph  has  been  shown.   It  also 
has  been  shown  that  computer  programs,  using  an  algebraic 
manipulation  language  FORMAC  can  be  used  in  the  subsequent 
evaluation  of  the  results  of  the  pattern  recognition  problem. 
It  is  possible  to  evaluate  in  symbolic  form  certain  system 
functions,  for  example  the  transfer  function  between  specified 
nodes  and  the  characteristic  equation  of  large  systems. 

A  FORMAC  program  can  also  be  used  to  advantage  in  a 
modified  root  locus  method. 

It  is  felt  that  the  application  of  an  algebraic  manipu- 
lation language  (FORMAC)  which  has  a  wide  range  of  capabilities 
as  compared  to  other  more  restricted  languages   which  may 
handle  only  polynomials,  is  of  great  value  in  the  solution  of 
algebraic  problems  in  control  theory  and  in  network  theory. 

Of  particular  importance  is  the  fact  that  FORMAC  is  a 
superset  of  PL/1.  All  the  capabilities  of  PL  1  are  usable 
together  with  FORMAC. 
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COMPUTER  PROGRAMS 

1.  Program  DTERM,  a  FORMAC  program 
Purpose:  To  evaluate  a  determinant  of  a  square  matrix 

of  order  up  to  9 .   The  elements  of  the  matrix 
are  univariate  polynomials  with  symbolic 
coefficients.   Given  is  one  example. 

2.  Program  FINDLOOP,  a  PL/1  program 
Purpose:  To  evaluate  a  control  vector  LIST1(4)  which  is 

needed  for  the  Mason  flow  graph  analysis.   Given 
is  one  example. 

3.  Program  MAS0N_1 ,  a  FORMAC  program 
Purpose:  To  find  sets  of  touching  loops  from  given  topo- 
logical matrices  describing  a  Mason  flow  graph 
and  to  evaluate  the  graph  determinant.   Given  is 
one  example . 

4 .  Program  MAS0N_2 ,  a  FORMAC  program 

Purpose :  To  evaluate  the  cof actors  in  the  Mason  gain  for- 
mula.  Given  is  one  example. 

5.  Program  RTLOC ,  a  FORMAC  program 

Purpose:  To  find  the  equation  of  the  root  locus  from  a 

given  characteristic  equation  in  symbolic  form. 
Given  are  three  examples. 
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RTLOC   Output,    Cycle    1 
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RTLOC  Output,  Example  1 


ro 


m 


c\j 


oo 

♦ 

i 
oo 

00 

+ 
I 

00 

+ 

+ 
OO 

o 

r-1 

II 
& 

lu 

I 

o 


o 

I 

< 

2. 
O 
t— i 
OO 


c\j 


<\i 
I 

< 

t— I 
00 

o 

LU 

2 
o 


ro 


<: 

¥ 
O 
►—I 

OO 

o 

+ 

II 

& 
o 


i  •— * 

1       h~   1 

i  % 

I    r-*    1 

i  ^ 

1     1   1 

i  + 

1  ^    1 

1    <  1 

1  00 

t  2.  1 

i   < 

1  0  1 

1   2 

1  *~~*  t 

1   o 

1   00  1 

•  »-h 

1   oo 

1   +  1 

1  Ol 

1  01   1 

1    < 

1   <  1 

1   O 

1   2:  1 

I    LU 

1   0 1 

1    21 

1   ►— <  1 

1   O 

1  00  1 

1  o 

00  1 

1 

+  1 

1  OJ 

ro   1 

1  < 

<  1 

1   21 

2:  1 

1   o 

O  1 

1   *~ ' 

•— •  1 

1   oo 

00  | 

1   < 

Is-  1 

1   o 

■— •  1 

1    LU 

1    2T 

+  1 

1    O 

|     *-H 

ro   1 

1      * 

<  » 

O  1 

1   -4  1 

LU  1 

1    OJ  | 

Si   | 

O  | 

1   +  1 

•— 1  j 

1   <  1 

*  t 

1   O  | 

1    LU  | 

00  1 

•   2:  1 

1   0  1 

•  • 

t   •— 1  1 

1   %  1 

ro   1 

<  1 

1   0 1 

2T  1 

j        T-i      | 

O  l 

»-H  | 

1        +   1 

00  1 

1       <■       1 

<■  \ 

1        5    1 

0  » 

1        ^   I 

LU  | 

1        ' '   1 

3E  t 

1        OO   | 

O  1 

<  I 
o  I 

LU  I 


o 


+  I  * 


r\j 


I 

<  • 
O  I 
LU  I 

21   t 

o  I 


Od 


ro 


(M 


< 
O 

O0 

+ 

< 

2! 
O 
1— 1 
00 


(\j 


< 
o 

LU 

2 
O 


< 

r 
o 

00 

ro 


O 
O 

_i 
cc 


o  I 

-•  I 
I 

+  I 
I 


86 


RTLOC  Output,  Example  2 
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RTLOC  Output,  Example  3 
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RTLOC  Output,  Example  3 
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RTLOC  Output,  Example  3 
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RTLOC  Output,  Example  3 
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